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Schur's transforms of a polynomial are used to count its roots in the unit disk. 
These are generalized them by introducing the sequence of symmetric sub-resultants 
C/3 ■ of two polynomials. Although they do have a determinantal definition, we show that 

they satisfy a structure theorem which allows us to compute them with a type of 
Euclidean division. As a consequence, a fast algorithm based on a dichotomic process 
and FFT is designed. 

We prove also that these symmetric sub-resultants have a deep link with Toeplitz 
matrices. Finally, we propose a new algorithm of inversion for such matrices. It has 
the same cost as those already known, however it is fraction-free and consequently 
CN \ well adapted to computer algebra. 



1 Introduction 



Let P = a + a\X H h a d X d be a polynomial in C[X]. In 1918 Schur gave 

a method to compute the number of roots of P in the unit disk [28]. This work 
was completed by Cohn in 1922 [7]. 

The so-called Schur-Cohn algorithm works as follows. Suppose that ciodd ^ 
and define the reciprocal of P by P* = X d P(l/X). Compute the following 
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sequence of polynomials : 

T(P) = P(0)P - lc {P)P*, T k {P) = T{T k ~\P)), 



where lc (P) denotes the leading coefficient of P. This sequence is finite : it 
has at most deg(P) polynomials with decreasing degrees and real constant 
terms. It is the variation of the signs of these constant terms, all supposed to 
be non-zero, which gives us the number of roots of P in the unit disk. See 
Henrici [17] or Marden [25] for a precise description of this algorithm. 

In this primary version, two difficulties arise. First, the algorithm does not 
work for every polynomial : if the difference of the degrees of two successive 
transforms T k (P) is more than one, or if some constant terms are zero, it is 
not possible to compute the number of roots of P. Second, the exact compu- 
tation of these transforms suffer from an exponential increase of the size of 
the coefficients : at each step, the length of the coefficients is approximately 
doubled. 

For these two reasons, we introduced the new sequence of Schur-Cohn 
SUBTRANSFORMS (see Saux Picart [33]). These subtransforms are equal 
to T k (P) up to a multiplicative factor, can be computed for every polynomial, 
have a determinantal definition, and an approximately linear increase is their 
coefficients. Moreover from the constant terms, we can compute the number 
of roots of the polynomial in the unit disk, using an adapted rule of signs. 

Later on, it appeared that the sequence of the Schur-Cohn subtrans- 
forms is linked to the sequence of the successive remainders of P and P* in 
a special "symmetric" division (see Brunie and Saux Picart [5]). This di- 
vision consists in eliminating from the largest polynomial as many monomials 
as possible from the top as well as from the tail by adding good multiples of 
the "divisor" . In the article cited above, we give a structural theorem, which 
describes the link between these two sequences built from P. 

In the present article we generalise the definition of the Schur-Cohn sub- 
transforms and the symmetric division of two polynomials to a general 
situation (no restriction on P and P*). We will speak of symmetric subre- 
SULTANTS of two polynomials. We are then able to formulate a new general 
"structure-theorem" which constitutes a central result of our work. With this, 
we compute the sequence of symmetric subresultants, using a EuCLlD-like 
algorithm instead of the determinantal definition. A dichotomic process and 
DFT allow us to produce a fast algorithm. Our methods are adaptated from 
ideas introduced by Schonage for the computation of Euclidean remainder 
sequences in [29], and by Lickteig and Roy in [23] for the computation of 
classical subresultants. The algorithm cost is of 0(M.(d) log d) arithmetical 
operations, where M. (d) denotes the cost of the multiplication of two polyno- 
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mials of degree d. 



We will not describe the application to the number of roots of a polynomial 
in the unit disk as it has already been discussed in [5]. However there are 
well-known relations between the problem of root isolation and TOEPLITZ 
matrices (see for example, M.G. Krein and M.A. Naimark [20]). We use 
these links to give, in the last part, a fast algorithm for solving Toeplitz 
systems with exact computation. It has the same cost as the well-known al- 
gorithm of Brent, Gustavson and Yun in [2], or those of Gemigniani in 
[13]. Morover, it is fraction free and consequently well adapted to computer 
algebra. We also give a new way to compute the signature of a Hermitian 
Toeplitz matrix. 

This paper is organised as follows. Section 2 introduces notations and defini- 
tions. In Section 3, we state the structure-theorem. Section 4 describes how to 
efficiently compute the symmetric subresultants and the last section applies 
these results to Toeplitz matrices. 

Finally, we wish to thank M.-F. Roy and T. Lickteig for their help and 
interest in this work. 



2 Definitions and Notations 



Consider a subring D of C and define the valuation of a nonzero polynomial 
P E H>[X], denoted by v(P), as the greatest integer v such that X'° divides P 
(it is also named "X-adic valuation" in many books). For the zero polynomial 
put deg(0) = — oo and v(0) = oo. Denote by W the quotient field of D. 

We write co&(P) for the coefficient of order k of P. If degP = d, the leading 
coefficient co^(P) is lc (P) and the trailing coefficient co„(p)(P) is denoted by 
tc (P). Remark : if v(P) ^ 0, tc (P) is different from P(0). 

We will use Euclidean division of a polynomial A by a polynomial B in U>[X] : 
the notation quo(A, B) stands for the quotient and rem (A, B) for the remain- 
der; they have their coefficients in the fraction-field W. We say that the division 
is exact if quo(A, B) and rem(A, B) are elements of D. Please note : our defini- 
tion of exact division differs from another definition common in the literature 
where exact division simply means vanishing of the Euclidean remainder. 

Now, let us introduce the main object of our article. 
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2.1 Symmetric Subresultants 



Let A = Eto a i xi and B = T,i=ob i X i be tw0 polynomials in B[X]. We 
suppose that one of them at least, say A, has its degree equal to d ; B can 
also be formally considered as having degree d : if deg B = d! < d, B will be 
replaced by 0X d + ■■■ + 0X d - d ' +1 + B. We also assume that the valuation is 
for at least one of them, otherwise we divide both polynomials by a power 
of X to ensure this condition. Define : 



Sylv ? .(A£) = 



/ ct 



a 

b n 



bo 



ad 



bj 



d+j 



to be a submatrix of the full Sylvester matrix Sylv d (A, B). 

For £ — 0, . . . , d — j, let Syfv^ = Sylv^(A, B) be the following 2j x 2j square 

submatrix of Sy\Vj(A, B) : 



/a 



Sylv 



a,j-2 Uj-i+e a d 



do <2e + i dd-j+2 

a i dd-j+i 
bj-2 bj-i+e bd 



bo 



b. 



bd-j+2 
bd-j+i 



ad 

■ ■ a d -i a d 



b d 

bd-i b d ) 



} J 



3-1 



The sequence (<Sj)-i<j<d of symmetric subresultants of A and B is defined 
by : 



• = A, 

• S = B, 



Sj(A, B) = Ytl det(Sylv^)^, if 1 < j < d. 



Clearly, Sj is an element of 3[X] for any j. The last one, Sd is just the resultant 
of A and B. In the generic situation, Sj is of degree d — j and valuation 0. 
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However, the real degree could be less than d—j and the valuation greater than 
0. In order to describe these situations, we introduce the following definition. 

Let (a,P) be such that : 



v{Sj 







deg(Sj) = d-j 



and < 



v(S j+1 ) 
deg(S j+1 ) 



a 



d-j-0 



we will then say that the pair (Sj, Sj+i) is (a, /3)-defective. The case (0, 1) is 
just the general situation without special deflation. 

Just as for the classical subresultants, we can express the Sj through a Bezout 
relation between A and B. This is estblished in the next lemma. 

Lemma 1 Let A and B be two polynomials in U>[X] of the same degree d and 
valuation 0. For every j e {0, 1, . . . , d — 1}, there exist two elements in Tt>[X], 
Uj and Vj, such that : 

X j S j+1 = UjA + VjB. 

The degrees of Uj and Vj are at most j . These polynomials are unique under 
such an assumption. 

Proof : Using a matrix with coefficients in ©[X], we can write X^Sj+i as a 
determinant in the following way : 



X j S 



a 



b ■ 



dj-i X j aj + ... + X^a^ 
\ Xidj-i + ...+ X d - X a d _ 2 



a Xia x + . . . + X*- 1 * 



ad 



Xiao + .-. + X*- 1 - 



a d-j+l '■ 

l a d _j ad-j+i 

bj-! Xlbj + .-. + X^bd-! b d 

: x'Vi + • • • + -V" h,i i i 



... a d 



b X% + . . . + x d -%_ 



j'+i 



X% + ...+ X^ba-j b d _ j+1 ...b d 



We do not change the value of this determinant by adding to the (j + l)-th 
column a linear combination of the other ones. More precisely, call Cj the i-th 
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column of the above matrix (i = 1, . . . , 2j + 2). Then add to the (j + l)-th 
column d + XC 2 + ...+ X^Cj + X d C j+2 + ... + X d+3 C 2j+2 . We obtain : 



X 3 S; 



do • • • flj— i ^4 

••. : ; 

a Xi~ x A \ 
X*A a d -j 



■ ■ a d 



b ... bj-i 



B 
XB 



b X'-^B : 
X>B 



Expand this determinant according to the (j + l)-th column, putting A as a 
factor in the first j + 1 lines and B in the last j + 1 : therefore there exist two 
polynomials, Uj and Vj, of degree at most j such that : 

X j S j+1 = UjA + VjB. 

Furthermore, we can express these polynomials as determinants. We have : 



and : 



CLq . . . CLj — i 1 (Xd 

■•• : X ': 



a Xi- 1 \ 

X 3 ad-j . . . ad 
bo ... bj-i b d 



bn 



b, 



d-j ■■■ "d 
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a ... a 7 _i a d 



a : 

a d _j 

b ... bj-i 1 b d 

■•. : x ; 



b Xi- 1 \ 

X 3 b d -j . .. b d 



□ 



For j = 0, we simply have Si = b d A — a d B, i.e. U = b d and V = —a d . Using 
the determinantal definition of Uj and Vj, we see that : 

UjiO) = b ■ co^-OS,-), cojiUj) = lc (Uj) = b d ■ 5,(0), 

and also : 

VjiO) = -oo • co^(S,), oo,-^-) = lc (Vj) = -a d ■ Sj(0). 

Finally, we can observe that these polynomials are uniquely determined, first 
when A and B are co-prime, and then in the general case. (The proof uses the 
same arguments as for the extended Euclidean algorithm for polynomials; see 

[11]-) 



2.2 Symmetric division of polynomials 



The division we use is justified by the following lemma. 

Lemma 2 Let A, B e B[X], with B ^ 0, deg A = d > deg B = d - (3 and 
v(B) = a. There exist Q,R e D'[X] 7 where W is the fraction field of D, 
uniquely determined, such that deg Q = a + (5 and deg R < d — (a + (5), and : 

A = Q§- a + X?R. 

Proof : We sketch how to compute Q and R. First divide A by B/X a with 
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increasing powers of X up to order f3. We obtain : 

A = Q 1 ^- + Xf i R u 

with degQi < (3 and degi?i = d — (3. Then, compute the Euclidean division 
of Ri by B/X a : 

Ri = ^2^- + R, 

where deg R < d — (3 — a and deg Q 2 = a. Then, define Q by Q = Qi + X l3 Q 2 
to establish the claim. Uniqueness is proven as usual. □ 

The polynomial Q is called the symmetric quotient of A by B, noted squo (A, B) 
and R the symmetric remainder, denoted srem(A, B). 

It is clear that the computation of such a division has the same arithmeti- 
cal cost as ordinary Euclidean division. It requires, at most, d(a + (3 + 1) 
arithmetical operations. 

Historical note : We can find various kinds of "symmetric" division introduced 
by authors with specific aims. See for exemple, Jezek [19], Demeure and 
Mullis [9] . However, our definition is different from the one in [19] and, when 
a — (3, coincides with the one given by DEMEURE and MULLIS only in the 
case. 



3 Structure-Theorem for symmetric subresultants 

We now describe the relationship between the sequence of symmetric subre- 
sultants and the sequence of symmetric remainders of two polynomials. Our 
main result is : 

Theorem 3 Let D be a subring of C, and let A and B be elements of"D[X] 
of degree d and valuation 0. Let (Si) <i<d be the sequence of symmetric subre- 
sultants of A and B. Suppose that the pair (Sj,Sj +1 ) is (a, (3)- defective. We 
have : 

(1) • if a > and (3 > 1, then Sj + k = for k = 2, . . . ,a + (3 — 1 
• if a = and (3 > 1, then, if j > : 

5,(0) • S j+k = ^(Of-^+i for k = 2, . . . , (3 - 1. 
If j = 0,S k = S 1 (0) k - 1 S 1 for k = 2, . . . , (3 - 1. 
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• if a > and (3 = 1, then if j > 1 : 



lc (Sj)"- 1 ■ S j+k = (-l) fc lc (Vi) fe_1 • |j£ for fc = 2, . . . , a. 



If j = 0,b k d -S k = (-l) fc • lc (S 1 ) k - 1 . . . S 1 X- k+1 for k = 2, ■ ■ ■ , a. 
(2) In all cases, if j > 0, we have : 



lc (s.r ■ sm"- 1 ■ s j+a+ p = (-i)^ • lc (s J+1 r • tc (s^f- 1 ■ 

and if j = 0, then : 

hi . s a+ ? = . 6g . ic (s 1 r • tc (s^- 1 ■ 

(3) In all cases, if j > ; we have : 



lc (Sj) ■ Sj(0) ■ S j+a+(3+1 = -lc (S j+1 ) ■ S j+a+/3 (0) ■ srem(S i; S j+1 ) 

= -srem (lc (S j+1 ) ■ S j+a+ p(0) ■ Sj, S j+1 ) 

and if j = then : 

hd ■ S a+ p +1 = -srem (lc (Si) • S a+/3 (0) ■ S , Si) . 

One remarkable fact is that the last symmetric divisions are exact in D, as we 
shall prove later. 

Observe that Si can also be expressed as a symmetric remainder : by Lemma 
I, we have Si = b d A — a d B = srem(S_i, So). 

It could be helpful to the reader to visualize the different situations. 
(I) Case (Sj,Sj+i) defective on "each side", a > 0,(3 > 1 : 
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Sj-i 
Sj+i 



Nullity 

Sj+a+13 

Sj+a+(3+l 



(2) Case (Sj, Sj +i ) defective on the "right-hand side", a — 0, (5 > 1 : 



Sj 

Sj+i 



Sj+a+13 
Sj+a+f3+l 



D — Proportionality 



(3) Case (Sj, Sj+i) defective on the "left-hand side", a > 0,(3 — 1 : 
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DfX] — Proportionality 



S 



s 



'j+a+/3+l 



Proof : Roughly speaking, we can say that the rows of Sylvi(A, B) are made 
of A,XA,..., X l ~ l A, and B,XB, ..^X^B, identifying the vectors of the co- 
efficients of these polynomials with the polynomials themselves. Furthermore, 
we consider them all of formal degree d + % — 1. 

Preliminary work : By Lemma 1, we know the existence of two polynomials, 



As the pair (Sj, Sj + i) is (a, (3) -defective, Sj(0) and cOd-j(Sj) = lc(Sj) are 
different from zero. Because of the determinantal definition of Uj (see proof 
of Lemma 1 ) , we have : 

- if j > 0, m = b ■ lc (Sj) + 0, and Uj = b d ■ Sj(0) ^ 0, 

- if j = 0, m = Uj = b d ^0. 
Then, for every £ > 0, we have : 



Uj = Yji=oUiX % and Vj = 521=0 ViX % , such that : 



X 3 S j+1 = UjA + VjB 



j j 



= Y J U n {AX n ) + Y,Vn{BX n ). 



n=0 n=0 



J 



3 



X 3+e S. 




u n AX n+E +Y,VnBX 



(t) 



n=0 



n=0 



with m and Uj different from 0. 
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For k > 2, and i fixed between and k — 1, we can replace the (i + l)-th 
row of Sylv J+fc , XM by the linear combination of the rows X i A,...,X : > +l A 
and X l B, X'j +l B described in (t). For i = i we obtain X^ +l Sj + i on the 
(i + l)-th row of Sylv J+fc instead of X % A. The minors of order 2(j + /c) of 
this new matrix are equal to u times the corresponding ones in Sylv J+A; . This 
operation will be called the (i, ^-transformation of Sylv J+fc . The downward 
arrow means that the j rows directly below the {i + l)-st row are used. 

We define also the (j + i, ^-transformation for i = 0,...,k — 1 : this re- 
places the (j + % + l)-st row by X j+t Sj +1 which is a linear combination of the 
rows X'A, X j+i A and X*B, ...,X j+i B, by (f). In this case the values of the 
minors of order 2(j + k) of Sylv J+jfc are multiplied by Uj. 

We use these two transformations in four different situations, described below. 
For each, we have drawn the corresponding matrix resulting from Sj + k- on 
each diagram, the rows with large dash patterns delimit the j + k — 1 first 
columns and the j + k last ones needed for the computation of Sylv^ e (£ — 
0, d — j — k). The shadowed triangles highlight the coefficients of the matrix 
needed for the computation of Sj(0). 

We consider now the four different cases. 

• 1 < P < a. Two situations have to be distinguished. 
O If 2 < k < a, we use k (i, |)-transformations for i = 0, k — 1 in this 
order. We obtain the matrix Mx (fig. 1). 



k 



j + k 



a-k + 1 






^\ 1 i 










A 





j+k -I d-j-k+1 j+k 



Fig. 1. Shape of the matrix Mi 

For each £ e {0, ...,d — j — k}, the minor det(Sylv J+fc ^) of Sylv J+jt is 
equal to the corresponding minor of the above matrix divided by Uq. If we 
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denote this minor by dj +ky £, we have 



Uo ■ det(Sylv j+M ) = d j+k/ . 

O If a < k < a + (3, we use a (i, j ^transformations for i — 0, a — 1, in 
this order, and then k — a (j + i, ^-transformations for i = k — 1, a, 
again in this order. We obtain the matrix M 2 (fig. 2). 



j + k 




j+k-l d-j-k+1 j+k 

Fig. 2. Shape of the matrix M2 

With the same notation as in the first case, we have : 
< • u k ~ a ■ det(Sylv J - +M ) = d j+k/ . 
» < a < p. Once again two situations occur. 

O If 2 < k < j3, we perform k (j + i, f ^transformations with i = k — 1, 0, 
in this order. We get the matrix M 3 (fig. 3), and we have for £ e {0, d — 
j-k}: 

u) ■ det(Sylv j+M ) = d j+k/ . 

O If /3 < k < a+(3, we use (3 (j+i, |)-transformations with i = k — 1, k—f3 
in this order, and k — (5 (i, |)-transformations with i — 0, k — (3 — 1, in 
this order. We get the matrix M 4 (fig. 4), and : 

u k ^ ■ u] ■ det(Sylv i+fc/ ) = d j+k/ . 



We now prove the theorem, step by step. 
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\ 


min(fc, a + 1) J 


^ 


^^^^ 



j + k-l d — j — k + 1 j + k 

Fig. 3. Shape of the matrix M3 



N \ (3 +|cT- k 






a+lj 




2(3 -k 








j + k-l d-j-k + 1 j + k 

Fig. 4. Shape of the matrix M4 
Proof of (1) : 2<A;<a + / 3-l 



Since we have to show the nullity of Sj + k for k = 2, a + (3 — 1, we need to 
show that the coefficients det(Sylv j+k,e) vanish for £ — 0, d — j — k. This is 
equivalent to showing that d j+k ^ = 0, for one of the matrices M l5 M 2 , M 3 or 
M 4 , because uo and Uj are both different from zero. 

• Case a > 0,(3 > 1 
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Suppose that 1 < (3 < a and 1 < k < a. We use M x : the submatrix 
corresponding to dj +kt g has at most one nonzero element on its first row. We 
use the corresponding column to expand it. The first row of the remaining 
minor has only zeros since (3 > 2. Hence dj +k / = 0. 

If 1 < (3 < a < k < a + (3 - 1, we use M 2 . We have a + (3 - k > 1 and 
then 



[(2a-k + l)+0\ -a>2. 

It follows that there are at least two among the first a rows for which 
at most one entry is nonzero, namely on the (j + k)-th column. Developing 
dj +k j along those two rows shows that it is zero. 

If 1 < a < (3 and 1 < k < (3, we use M 3 to expand dj +k ,e along the (j+k)- 
th row, which has at most one nonzero coefficient. As min(A;, a + 1) > 2, 
the row immediately above also has this property, and we get dj +k / = 0. 

Finally, ifl<a<j3<k<a + j3 — 1, we use M 4 . Once again, in 
dj+k,e we have two successive rows with only one non-zero coefficient, on the 
(j + k)-th column (because (a + 1) + (2(3 - k) > (3 + 2). 

In every case, we see that, if a > and f3 > 1, then Sj +k = 0. This 
establishes the first part of 1. 
Case a — 0, (3 > 1 

As 2 < k < a + (3 — 1, we have 1 < k < (3, min(/c, a + 1) = 1 and we can 
expand the minor dj +k ,e, using the rows j + k, ...,j + 1 in M 3 , in this order, 
and then, using the last k columns. We obtain, for every £ — 0, d — j — k : 

d j+k4 = co e (S j+1 ) ■ tc (^■ +1 ) fe - 1 • b k d ■ 5,(0). 

(The factors are written from left to right, in their order of appearance in 
the successive expansions.) As dj +kt g = u k det(Sylvj +k /) and Uj = bdSj(0), 
we have : 

Sj(0) k 1 ■ Sj +k = tc (Sj + i) k 1 • Sj+i. 
If j = 0, Uj = b d and Sj(0) does not appear in dj +k /- Hence : 

5 fc = tc(S' 1 ) fc - 1 .S' 1 . 

a > 0,(3 = 1 

We have 2 < k < a and we use M 1; expanded along the first k columns, 
and then along the first k rows. We obtain (the factors appear in order of 
expansions from the right-hand side of the formula) : 

d 3+k/ =(-l) k ^ +k+2) ■ b k ■ (-l) k « +k+ %o J+k _ 1+e (XiS 1+1 ) 

= (-l) k ■ b k ■ co k _ 1+e (S j+1 ) ■ lc (^■ +1 ) fc - 1 • lc fa). 

The result follows. If j = 0, the computation is the same : however, in 
this case, all the rows of block B collapse. 
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Proof of (2) : k = a + (3 



If (a, (3) = (0, 1), the result is trivial. So, we suppose that (a, (3) ^ (0, 1). For 
j 7^ 0, we distinguish two cases. 

• (3 < a 

We use the matrix M 2 and expand it along the row of order (3 to obtain 



d j+k ,e = < " u j ' det(Sylv j+k , e ) 



4 



<-^-(-ir-co a+ ,(^ + i)-A, 



where no = j + a and A is a minor independent of £. 

Then, we expand A along the first (3 — 1 rows, and see that : 

A = (-ir-tc(S', +1 )' 3 - 1 -A 1 , 

with n-i — (j + a){(3 — 1). We continue expanding A x along the first a — (3 
rows ; we have : 

A 1 = (-ir-lc(S' J+1 ) a -' 3 -A 2 , 

with n 2 = (j + a)(a — (3). We can then use rows j + 1, j + (3 to compute 
A 2 : 

A 2 = (-ir-lc(S' J+1 )' 3 -A3 

(n3 = af3). Finally, A 3 can be expanded using the first a columns and the 
last (3 ones : 

A3 = (-ir-&o-^(o), 

with 114 = ja. In summary, we have obtained : 



u -Uj 



det(Sylv j+k , e ) = (-if -b^-tc (Sj^-lc {S j+1 ) a 'S j (0)-co ct¥t (S j+1 ), 



with N = n + rii + n 2 + n 3 + n 4 = a(a + (3) mod 2. As this computation 
is valid for every £ — 0, d — j — k, we have : 

lc (Sjf ■ SM* 3 - 1 ■ S j+a+P = • lc (S i+1 ) a • tc {S J+1 f- 1 ■ 

a < (3 

We use the same method as in the previous situation, starting with M4. 
We expand it along the row of order j + (3 and obtain : 



dj+kjt = < • u j ■ det(Sylv w ) = u% ■ u] ■ (-l) n ° ■ co a+e (S j+1 ) • A'. 
We expand A' along its rows j + (3 + 1, + k to obtain : 
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then again along its rows j + (5 — 1, j + a + 1 to obtain : 

AiH-ir-tc^y 3 -"- 1 -^, 

and then along its first a rows : 

A' 2 = (-irs. t c(Vi) Q -A^ 
to finally find that A 3 = A3. We now have : 

n' = a, 

1 2 
n 1 = a , 

n' 2 = a(/3 — a — 1), 
n 3 = a U + 

We obtain exactly the same final relation as in the case (3 < a. 

If j = 0, we have uq = Uj = bd, and Sj(0) disappears at the end of the 
successive expansions of the minors. Therefore we get : 

b a d ■ S a+ , = (-!)«(«+« -fiff-Ic (S^ • tc (S'i) /3_1 • 



Proof of (3) : 

Here we cannot use the same transformations of Sylvj +a+ /3 + i as above. 

We suppose first that j > 0. Let R = — srem(5j, Sj+i) and Q = squo (Sj, Sj+i). 
There exist four polynomials Uj-i, Vj-i, Uj and Vj such that : 



X j - 1 S j = Uj- 1 A + Vj- 1 B, 
X j S j+1 = UjA + VjB 

with deg(^-_i) <j-l, deg(y j _ 1 ) < j - 1, deg(^) = deg(V-) = j. We also 
have : 

and deduce that : 

Xi+ a +PR= (QUj - X a+1 U^ X )A + (QVj - X^V^B 
= UA + VB. 
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As deg(QL r i ) = j+a+(3 and deg(X a+1 Uj-i) < j+a, we have deg U = j+a+(3. 
Likewise, deg V = j + a + (3. 



Also : 



\c(U) = lc (Q)lc(Uj) 



lc^) •6 d -g J -(0) 
lc 



The equation X j+a+/3 R = UA + VB, with deg(C/) = deg(V) = j + a + shows 
that X^ +a+l3 R can be obtained as a linear combination of rows of Sylvj +a+ p +1 . 
As in the previous steps, we transform the row j + a + f3 + 1 and obtain a 
matrix which has the following structure : 



^-X3+ a +f 3 R 

^Xj+c+l 3 B 



Therefore, for £ — 0, 1, d — (j + a + f3 + 1), we obtain : 

a ■ ■ ■ a J - +a+/3 _i aj +a+ p + e a d 



lc(C/)-det(Sylv j+Q+/S+1 ^) 



b 



ao a t~\ a d-j-a-/3+l 

co e (R) 

bj +a +i3-i b j+a+/3+ e b d 



■ a d 

■ 



bo 




b d - 



■j-a-P 



j+»+0 



Expanding these determinants along the last column, and then along row 
(j + a + (3 + 1), we see that : 

lc (U) ■ Sj +a+/ 3 + i = b d • Sj +a+/ 3(0)R. 



We use the value of lc (U) already computed to obtain the desired result : 



lc (Sj) ■ Sj(0) ■ S j+a+/3+ i = lc (S j+1 ) ■ S j+a+ p(0)R. 



b d 



When j = 0, the polynomials Uj-i, Uj and Vj are very simple, as we 
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have : 



S = 1.B, S 1 = b d A-a d B. 



The expression of lc (U) is now : lc (U) = 
computation is unchanged, and we obtain : 




. However, the rest of the 



lc (S )-S a +13+1 — 



□ 



Remark : If we define the Toeplitz-Bezoutian of two monic polynomials 
P and Q of the same degree as the matrix Bez(P, Q) whose entries are the 
coefficients of the polynomial 



If sCi(M) denotes the i-th ScHUR-comp lenient of the square matrix M when- 
ever it exists, one can see that we have : 



(See Bini and Pan [3] p. 169 for the classical result over the Euclidean remain- 
der sequence. Proof uses same methods). 



4 Computation of the Symmetric Subresultants Sequence 

The previous theorem gives us a direct method to compute the sequence of 
symmetric subresultants of two polynomials A and B, of same degree d and 
same valuation 0. It uses symmetric divisions instead of the determinantal 
definition. With parts 2 and 3 of Theorem 3, we can compute the subsequence 
(5 , fcJi=o,...,s (s < d) of the sequence of the symmetric subresultants, such that, 
for each index i, the pair (S^, S^+i) is (ctj, /^-defective. This implies that, 
for each i, S ki is of valuation and degree d — ki (we have k = as So = 
B). Denote by Qi the i-th. symmetric quotient of (S^, 5^+0 ■ The sequence 
(Ski)i=o,..., s is obtained by the following Euclidean-like algorithm : 



P(X)Q*(Y) - P*(Y)Q(X) 
1-XY 



Si(0)\c {S i )sc i {Bez(S- 1 , S )) = Bez(Si, S i+1 ). 




X^\c(S kl )-S kl (0)-S, 



fc2 + l> 
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1c (S ks +i) ■ S ks+1 (0) ■ S ka = Q s 




For such an algorithm, a classical analysis of cost gives a bound of 0(d 2 ) 
arithmetical operations. In the important case of Z, we use HADAMARD's 
bound for a determinant : if the size of all the coefficients of the polynomials 
is bounded by a, then the size of the coefficients of all the S ki is bounded by 
2d(a + log(cT)). Therefore, in the case of Z, the binary cost of the algorithm is 
in 0(d 2 M(2d(a + logcf))) where M(t) denotes the cost of the multiplication 
of two integers of absolute value less than 2*. 

However, this algorithm can be improved. In a previous article (see [5]), we 
studied the case where B is the reciprocal polynomial of A. In fact the im- 
provement we gave can be applied to every pair of polynomials A and B in 
3[X] of same degree d and valuation zero. The next section is devoted to 
showing this. 

The ideas we develop here are adaptations to the case of symmetric subre- 
sultants, of ideas already known for ordinary subresultants (see [21], [22], [23], 



4-1 Transition Matrices 

One idea is to express the transition from a pair 
(S ki , Sfci+i) to a pair (S ki+1 , S ki+1 +i) with an appropriate matrix. 

Let A and B be two polynomials in D[X] of same degree d and valuation 0. 
Suppose the pair (Sj, Sj + i) to be (a, /9)-defective ; set k = j + a+j3 and denote 
by Q the symmetric quotient of lc (Sj + i)Sk(0)Sj by Sj + \. With formulae 2 and 
3 of the Structure-Theorem Th. 3, we can write, for j > : 



[27]). 




with 







( a+/3 ) a lc(5 J+1 )"tc(S J+1 )/ ; '- 1 \ 
lc(S' j )«S' i (0)' 3 - 1 ^ 



lc(Sj + i)S fc (0) ya+1 
\ lc(S,)5,(0) ^ 



(1) 



lc(S^(0) / 



In the case j = 0, we obtain : 
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with 



0,fc 







»->o 1 




= M 0tk ■ 


^X k S k+1 J 





^ ^fcg fcglc (Si)°tc (Si)/ 3 ' 1 ^g-l > 



lc(5i)g fc (0) Y « 
6 d 



(2) 



Furthermore, we have (for j 

( 




We can now state a general definition. 

Definition 4 Lei A = Z)f =0 a^X 1 and £? = Sf =0 ^ e ^0 polynomials of 
3[X] of same degree d and same valuation 0. Let (Si)^i<i<d be the sequence 
of the symmetric sub-resultants of A and B. We denote by (/cj)j=o,..., s (with 
k = < ki < ... < k s ) the sequence of indices such that (S ki ,S ki+ i) is 
(ojj, Pi) -defective. 

Then, for i,j G {0, s}, with i < j , we denote by M k . ;k . the matrix defined 
by : 



Mi 



Mi 



■Mi 



■j-2,kj-i 



■Mi 



where the matrices M kl)ke+1 are defined by the above formulae (1) and (2). If 
i > 0, we have : 



( x k >- l s k . 

\^X k iS kj+ i 



M k . fc . • 



1 ' x k ^s ki 



X ki S, 



ki+1 



and if i = 



x k ^s kj 

X kj S kj+ \ 




We call the matrix M kiik . the transition matrix from the pair (S ki , S ki+ i) to 
the pair (S kj , S kj+ i). We denote by M ki the transition matrix from (A,B) to 
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(S ki ,S ki+ i) : 



M ki = M , 




with the convention that M 0i o is the identity. 

We can now justify the assertion of the previous section : all the quotients 
(and remainders) involved in the Structure-Theorem are fraction-free. 

Proposition 5 Let A = Yn=o Q>iX l and B = Yn=o hX l be two polynomials of 
3[X] of same degree d, and same valuation 0. Let (Si)-i<j<d be the sequence 
of the symmetric sub-resultants of A and B. Let j e {1, d — 1} be such that 
(Sj, Sj+i) is (a, (5) -defective. Put k = j + a + (5. 

Then the symmetric quotient of lc (Sj + i)Sk{0)Sj by Sj+i belongs to V>[X], as 
does the symmetric remainder. 

Proof : By Lemma 1 we have for i > : 



X^S^U^A + V^B, 
X'Si+^UiA + ViB. 

Therefore, we obtain, for each j > 0, the following expression of Mj : 



Mi 



V Uj v 3 , 



We can directly deduce from (1) and (2) the value of det(M ?i fe). Moreover, if 
we consider the first line of 







f Xi- l S^ 




= M jik ■ 




^X h Sk+i J 


yx 3 s j+1 J 



we see that \c(S k ) = ( - 1 ) ^ a+ ^ a lc - ■ Therefore, we obtain, for 

j > : 

\c(S k )S k (O) xa+0 



and j = yields : 



det(M,- fc ) 



det(M 0lfc ) = 



lc(^(0) 
\c(S k )S k (0) xa+ ^ 
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As above, we denote by (/cj)o<j<m the indices such that (S ki , S ki+ i) is (c*i 
defective with k = and k m = j. We have : 



fm— i 



det(M j ) = det(M ) • II det(M fci)fci+1 ) , 



z=0 



m—1 



= -b d II det(M fcijfci+1 ), 



i=0 



6d M lc(^)^(0) 



m—1 



=-ic(^ m )5 fcm (o)x fei - fe »- 1 . n x^- fe % 



i=i 



-lc(5 fc .)5 fc .(0)X 



Consequently, the matrix Mj is invertible and we easily see that, if j > 



det(M i )Mr i = 



V* -v;--i 



When j = 0, we get : 6 rf M 1 



a d 1 
b d 



By definition of Mj, we have for < j < k, M j>k = M k ■ Mj 1 . Then for j > 
-]c(S j )S j (0)X*- 1 M j , k = 1 U 



£4-1 Vk-l 

u k v k 



v 3 -v^ 

-Uj ^-i 



and for j = : 



-b d M 0yk 



U k -i V k -i 

u k v k 



a d 1 
b d 



Identifying the bottom right-hand side entries of these matrices, yields 
if J > : 

X^Q = Uj-M - UkVj-! e B>[X], 
and Q = —U k when j = 0. □ 
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4-2 Symmetric truncation 



The computation of the symmetric quotient of two polynomials does not in- 
volve all of their coefficients. In fact, we only need the leading and trailing 
terms of the divisor. More generally the computation of successive symmetric 
quotients can be done with only the knowledge of a few leading and trailing 
terms of the first divisors. This way it appears cheaper to compute successive 
quotients instead of successive remainders, as we use only small parts, which 
we will refer to as "symmetric truncation" of the polynomials. 

First we define the symmetric truncation of a polynomial. 

Definition 6 Let P = EtoPi xi be an element of B[X], P ^ 0. For £ e 
{1, Ki/2j} ; we denote by P\e the polynomial 

P\e = Po + ---+ Pe-iX*' 1 + p d - £+1 X e + • • • + p d X 2i -\ 

For i = 0, we write P| = 0, and for £ > [d/2\ , P\ t = P. 

We now analyse the cases where truncation of two polynomials does not affect 
their symmetric quotient. 

Lemma 7 Let P and P\ be two polynomials o/B[X] such that deg(P) = d, 
deg(Pi) = d- (3<d, v(P) = and v{P 1 ) = a > 0. Then, 

squo (P, Pi) = squo (P|( Q +/3+i), Pi|( Q +/3+i)), 

where Pi is considered as a polynomial of degree d in order to compute its 
truncation. 

Proof : Set P = P|( Q+/3+ i), Pi = Pi|( Q+/3+1 ) and 7 = d — 2(a + (3) — 1. We 
have degP = 2(a + (3) + 1, degPi = 2a + (3+1, v(P) = 0, and v{P 1 ) = a. 
Then, let us consider the following symmetric divisions : 



P = q3- + X i3 R with deg(P) < d - a - (3, 

P = Q^- + X^R with deg(P) < a + P + 1. 
X a 



We have deg(Q) = deg(Q) = o. + (3 and we can write : Q = QiX 13 + Q 2 
and Q = Q\X® + Q 2 , where degQ 2 and degQ2 are strictly less than (3, and 
deg Qi = deg Q x = a. 

Then : 
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P-X^P = Q Pl x X a " Pl +(Q-Q)^ + X?(R - Xi R). 

Since deg(P — X 7 P) and deg(Pi — X^Pi) are less than d — a — (3, we see that 
deg(Q — Q)</3, and therefore, Q 1 = Q\. Similarly, we compare the valuation 
of both sides at the identity : 

P-P = Q^^ + (Q-Q)^; + X^R-R). 

As v(P-P) > a+f3 and v^-P^ / X a ) > (3, we conclude that v(Q-Q) > (3, 
i.e. Q 2 = Q 2 . Hence Q = Q. 

□ 

We can also compare the truncation of the symmetric subresultants of two 
polynomials with the symmetric subresultants of their truncations. 

Lemma 8 Let A and B be in B>[X] of same degree d and valuation 0. Let 
(Si)-i<i<d be the sequence of the symmetric subresultants of A and B. Let, 
(Sj)-i<j<2e-i be the sequence of the symmetric subresultants of A\t et B\e (£ 
fixed in {1, [d/2\}). Then for 1 < j < I, we have : 

S m-j) = S M-3)- 

Proof : The proof is based on the definition of the coefficients of the sym- 
metric subresultants. Set A = Y,i=o a iX l and A = Yh^o 1 a iX l ( respectively 
B = T,t=objX i and B = E^M"*)- For < k < £ — j, the coefficient of 
order k of is given by : 



a o ■ ■ ■ a j-2 a k+j-l a 2l-l 

d : : 
dk d 2 e-j ■ ■ ■ d 2 £-i 
b . . . bj_ 2 bk+j-i b 2 £-i 

b k b 2e _j . . . b 2 e-i 



c °k(Sj\e-j) = 
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QO • • ■ a j-2 O-k+j-l &d 

afc cid-j+i ■ ■ ■ ad 
b . . . bj_ 2 h+j-i b d 

b 

b k b d _ j+1 ... b d 



cOkiSjie-j). 



In the same way, if £ — j < k < 21 — 2j, we have : 



ao ■ ■ ■ a,k+2j-i a>2t-i 

a : : 

hk+j a 2 e-j ■ ■ ■ &2i-\ 

bo ■ ■ ■ bj-2 bk+j~i b 2 e-i 

k ! ! 

bk b 2 e-j ■ ■ ■ bu-\ 

ao ■ ■ ■ aj-i a d ~2t+k+2j a d 

a : : 

a d _ 2 e+k+j+i a d -j + i . . . a d 

bo ■ ■ ■ bj-i b d -2i+k+2j b d 

bo ■ ■ 

b d -2e+k+j+i b d -j+i ■ ■ ■ bd 
= cod-2e+k+j+i(Sj) = co k (Sj\(e-j))- 
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Therefore Sj\^j) and have the same coefficients. 



□ 



As a consequence, we have Sj\ k = Sj\ k for every k such that < k < I — j. 
Also Sj(0) = Sj(0) and lc (Sj) = lc (Sj) for every j < i. 

Further, for a given £, we can predict how many symmetric quotients will be 
preserved if we replace A and B by A\£ and B\£ in the computations. 

Theorem 9 Let A and B be in D[X] of same degree d > 4 and valuation 
0. Let (Si)-i<i<d be the sequence of the symmetric subresultants of A and 
B. For t G {2,..., [rf/2j} ; let (5j)-i<j<2£-i be the sequence of the symmetric 
subresultants of A\£ et B\£. 

Let (ki)o<i< s , respectively (ki)o<i< s >, be the indices such that the pairs (S ki , S ki +i), 
respectively (5~ , <% +1 ) ; are (aj, defective, respectively (Si, Pi) -defective (we 

have k = k = 0). 

For each i such that S ki +i ^ ; set Qi = lc (S , fc i+ i)S , fe i+1 (0)squo (S^, S ki +i) and 
for each i such that S^ +1 ^ 0, set Qi = lc (S^ . +1 )>% +1 (0) 
squo , Sj: Then M kuk . +1 , respectively M~~^, are the transition ma- 
trices of the sequence (5j)i<7<d; respectively (5j)i<j<2^-i- 

Let m be an index such that 1 < m < s and let k m + 1 < t, then for all 
i — 0, 1, m — 1, we have : 

C^i C\~i} Pi Piy Qi Qi} ^i ki, 

and finally, +i = M kiM+1 . 

Proof : First notice that for any i — 0, m — 1, we have k i+ i = ki + at + Pi, 
it follows that : 

m— 1 

J2®i + Pi<t- 

i=0 

For each j < i, by Lemma 8, we have Sj\£-j = Sj\£-j. Therefore, for each 
j = 1, 2, ...,£— 1, we have Sj(0) = Sj(0) as well as lc (Sj) = lc (Sj). Then, we 
see that ki = ki for every i = 0, 1, m. Furthermore, as ki + \ — ki = oti + Pi, 
and ki + i — ki = Sj + Pi, we have ctj + Pi = on + $ for every i — 0, 1, m — 1. 

We claim that ctj = Sj (i — 0, . . . , m — 1). This will also imply that Pi = Pi for 
each i = 0,1, ...,m- 1. Indeed, we have S ki+1 \£_ ki _i = S' fci+ i^_ fcj _i. Therefore, 
the £— fcj — 1 bottom coefficients of S^+i and S^+i are equal. But v(S ki+ i) = cti 
and we have ki + cti + Pi = k i+i < k m < t. Thus ctj is less than t — ki — Pi < 
t — ki — 1. The valuations of S ki +i and S^+i must then be equal. 

Having proved that the sequences of indices {ki) <i <m , {ai) <i <m , (Pi)o<i<m 
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are equal to their counterparts, we now show the equality of the symmetric 
quotients. 

First we have, by Lemma 7 : 

Qi = squo (lc (S ki+1 )S ki+1 (0) ■ S ki , S ki+1 ) 

= squo (lc(5' fei+ i)5' fei+1 (0) • S ki \ ai+/3i+1 , S ki+1 \ ai+ p i+ i), 

since (S ki ,S ki+i ) is (a iy /^-defective . 

If i < m and k m + 1 < £, we have ctj + fa + 1 < t — ki, and, by Lemma 8, 
S ki \ a . + p i+ i = S ki \ ai+ p i+1 . In respect of S ki+1 \ a . + p. +1 , the truncature is applied 
to S ki+ i considered of formal degree d — ki (Lemma 7). But, by Lemma 8, 
we have S^+i^+^+i = S ki+ \\ ai+ p i+ \, polynomials being truncated with their 
actual degree. However using formal degree d — ki instead of actual degree 
d — ki — Pi, we do not take into account so many coefficients and the equality 
of the truncatures holds as well. 

Since the leading coefficients and constant terms of the sequence (S^q^j^ and 
(Sj)o<j<£ are equal, we can write : 

Qi = squo (lc (S ki+1 )S ki+1 (0) • S ki \ ai+/3i+1 , 5^+11^+^+1) 
= Q l - 

Finally, inspecting the expression of the transition matrix M kuk . +1 given by 
(1) and (2), we see that all the ingredients have been proven to be equal for 
the two matrices M ki>ki+1 and M fciA+1 (i = 0, . . . , m - 1). 

□ 

4-3 Fast Algorithm 

We now describe the FSSR Algorithm which is written in pseudo-code further 
down. 

Let A and B be two polynomials of 3[X] of same degree and valuation 0. They 
are considered as global variables. The FSSR Algorithm takes as input a pair 
(S ki ,S ki+ i) of two successive symmetric subresultants of A and B, (a iy fa)- 
defective and an integer r < (d — ki). 

It returns the sequence of the symmetric quotients (Qj,ctj, j3j)i<j< v -i with v 
the largest index such that k v < ki + r. It returns also the transition matrix 
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In the general case, we are interested in finding the entire sequence of sym- 
metric quotients of A and B, and FSSR(S , , Si, d) with S = B, Si = 
lc (B)A — lc (A)B will suffice. This way, we compute the entire sequence of 
symmetric quotients except perhaps for the last one which can be obtained 
with an extra division. 

How does this work ? We use a strategy of divide and conquer, to compute a 
partial sequence at each step. Here is a description of each non-trivial step. 

Step 1 : If Ski+i is 0, we have already reached the end of the sequence of the 
symmetric subresultants of A and B. 

Step 2 : If r < 2, the algorithm performs symmetric divisions starting with 
the polynomials S ki \ r and S ki+ i\ r whose degrees are at most 3. It computes 
also directly the corresponding transition matrix. 

Step 4 : a call to FSSR (S ki \ r , S kt+ i\ r , is executed. 

Since the third recursive call, the coefficient of truncature is stricktly lower 
than L^y^J , and therefore Theorem 9 can be applied : the algorithm computes 
Qj, otj, Pj for j — i, . . . , u — 1 as well as M kitku , with u the largest index such 
that k u < h+ \r/2]. 

Step 5 : We compute S ku , and S ku+1 via M kijku . 

Step 6 : Then, via a symmetric quotient, we compute Q u and add it to the list 
of quotients already computed. M kuku+1 is computed as well as (S ku+1 , S ku+1+ i). 

This intermediary step is needed to guarantee that the coefficient of truncature 
in the next call to FSSR (step 7) is smaller than |~|] . 

Step 7 : We perform a second call to FSSR^S fcu+1 | r , S fcu+1+ i| r , r — (k u+l — kifj. 
We therefore obtain symmetric quotients Q u up to Q v -i with v the largest 
index such that k v + 1 < r + hi. 

Step 8 : We get together the pieces already computed. 

Remark : throughout the algorithm, instead of computing M kukm = M kjtkm ■ 
M kukj for < i < j < m < s, it is preferable to compute : 

M kt , km = (((lc (S kj )S kj (0)) ■ M k ^ km ) ■ M ki ^/(\c (S kj )S kj (0)) 

using the order of operations indicated by the parentheses. In doing so, we 
keep all computations in D[X] and the algorithm remains fraction-free. 
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ALGORITHM FSSR 



INPUT : - (S ki , Sfc 4 _|_i), a pair (on , /3; )-defective of symmetric 

subrcsultants of A, _B, 

- r a positive integer, r < d — ki. 

OUTPUT : - the list L := [Qj, a;, Pi, Q„_i, and M k . :kv , where v is the 

biggest integer such that k v < r + fej . 



MAIN PART : 1 - IF S fc . +1 = then RETURN L := [], and M := 7d 2 . 

2 — ELSE IF r < 2 then compute L using symmetric divisions of S ki \ r with 
S ki +i\ r an< i M kiku from definition. 

- ELSE 

3- r':= [fl; 

4 - Li :=FSSR(5 fc .| T .,5 fc . +1 | T .,r'); 

% Li contains : 

% Qj, Qj, ft, Qu_i, ct u -i, Pu-i, 
% we get also : M k .^ ku , 

% mt/i u, largest integer such that k u < r' + ki. 
5 - Compute S ku and 5 fcu+i by : 




M, 






- Qu = lc (S fcu+ i)Sfc u+1 (O)squo (5 fcu , S feu+ i) ; 

Li = L x \J{Qu}- M k ^ ku+1 = M kuiku+1 ■ M k . tkv 
Compute S ku+1 and S ka+1+1 by : 



7-L 2 :=FSSn(S ku+llr ,S ku+1+1 \ r ,r - (k u+1 - ki)); 
% L2 contains : 

% Qu+l,ct u +l,f3u+i, Qv-l, a„_i,/3t,_i; 
% u>e get also : M k l kv . 
% with v, largest integer such that k v < r + 
8 - L := Li (J L 2 ; M fc ., fc „ = M ku+likv ■ M k . iku+1 



END. 



We now consider its cost. 

Theorem 10 Let D be a sub-ring of C and let A and B be two polynomials 
of same degree d in D[X]. The algorithm FSSR(Sq, Si,d) with So = B and 
Si = lc (B)A — lc (A)B uses at most 

0{M(d)Aog(d)) = 0(dlog 2 (d) loglog(d)) 
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arithmetical operations in D (JA(d) denotes the cost in arithmetical operations 
of multiplying two polynomials of degree at most d in 3[X]). 

If A and B are elements ofZ[X] orZ[i][X], and if the size of their coefficients 
is bounded by a, then FSSR(S , Si,d) is executed in less than 

0((d 2 .(a + log(d)). log(d<7 + dlog(d)). log ( log(d<7 + dlog(d))). log(d)) 

binary operations on a multiband TURING machine, using DFT. 

Proof : Let us denote by CJ-{5) the cost in terms of arithmetical opera- 
tions of the computation of FSSR(S , Si, 5). We do not take into account 
the degrees of the polynomials S and Si, as, from the very beginning of the 
algorithm, these polynomials are truncated to order 5 and the degrees of the 
polynomials that we really manipulate are lower than 25—1. 

During the execution of FSSR(S , Si, 5), we use two calls of FSSR with 5 re- 
placed by |~|] . The intermediate computation consists of some multiplications 
and a symmetric division : the number of arithmetical operations is bounded 
by 0(M.(5)) . Therefore, we have : 



CF{5) < 2CT\ 



+ 0(M(S)). 



It follows that £F(5) is bounded by O (M(5) log(5)). Hence the first assertion 
with 5 = d. 



In the case of Z or Z[i], we follow the same arguments. However, we have to 
bound the size of the coefficients appearing in the algorithm. These coefficients 
are minors of Sy\v d (A, B). They can be bounded by Hadamard's formula : 
their size is less than r = 2d(o + log(rf)). The coefficients of the transition 
matrices M kitkj are of the same size. If A4(d,r) is the binary cost to compute 
the product of two polynomials of degree less than d with coefficients of size 
bounded by r, we get : 

CT{d,o) <0(M(d,T)\og(d)). 

This proves the result in the case of a multiband Turing machine. □ 

Remark : it might surprise the reader that we compute the sequence of sym- 
metric quotients instead of the symmetric sub-resultants. Indeed as far as 
applications are concerned the important elements are the symmetric remain- 
ders and not the symmetric quotients. In fact, the applications we know of 
use either the constant terms of a sequence of symmetric remainders, or a 
particular symmetric remainder. When the sequence of symmetric quotients 
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is known, the sequence of 5^(0) can be computed in O(d) as we can see in 
the introduction to Part 4. 

In this case, when a particular symmetric remainder is needed, computing the 
corresponding transition matrix is enough to determine this specific remainder, 
up to a few additional operations. 



5 Application to Toeplitz matrices 

In this section we consider the relationship between sequences of principal 
minors of a Toeplitz matrix and of the symmetric sub-resultants of polyno- 
mials. As a consequence, we will get new algorithms to compute the signature 
and the inverse of such a matrix. We do not improve the cost of algorithms 
presented in [2] and [13] and already used in the complex numerical case. 
However, in the case of integer coefficients, we control the size of results and 
use fraction-free computations; this is well suited for computer algebra. 



5.1 Relationship between Toeplitz matrices and symmetric sub-resultants 

We first establish a link between constant terms of the symmetric subresultants 
and principal minors of a TOEPLITZ matrix. 

Proposition 11 Let F = Y^i=ofiX l and G = Yn=o9iX l be two polynomials 
of equal valuation; we suppose that the degree of F is exactly d; the degree of 
G is formally considered equal to d but could be less. Let 

b i>l 

be the expansion around zero of G/F, and 

its expansion around infinity. Let Tk(F,G) = {U,j)i<i,j,<k be the Toeplitz 
matrix : 

Uj = Vj_i if i<j 
< t i:j = if i> j ■ 

U,j = u + v if % = j 
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Then, if (<Sj)-i<j<d £/ie sequence of symmetric sub-resultants computed with 
S-i = F and S = G, we have, for any k — 1, d : 



S k (0) = (-l) k .f*.f*det(T k (F,G)). 



Proof : As we have G = (—u — J2i>o u iX l )F, the following sequence of 
relations holds : 



90 = 

91 = ■ 



-uf - 

-uf\ - Uif 2 



- tirf-i/d-i - u d f d , 

- U d -if d , 



9d = -ufd 



Now define for k = 1, d, the following three k x k matrices : 



/ fd 0^ 

fd-l fd 



\ fd-k+1 fd J 



/ 9d o o^ 

9d-i 9d 



\9d-k+l 



9dj 



u, = 



u 

U\ u 



\Uk-l 



u 



Our relations can be translated by the following matricial relation : 



Gfc — —Uk • F fc . 



Likewise, comparing the coefficients of G — (v + J2i>o v iX l )F, we obtain : 



with 
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fo fk-l 
fo fk-2 



fo j 



V* 



9o 



9o 



Vk-2 



9k-i 

9k-2 

9o J 



These relations imply : 
( 




F k \ _ lF k F k 
F k I \G k G k 



(lit denotes the identity matrix of order k.) Now, we can compute the determi- 
nant of each side. For the left most matrix we subtract the i-th column from 
the (i + k)-th one (i = 1, k). The result follows. □ 



5.2 Signature of an Hermitian Toeplitz matrix 



Given an Hermitian TOEPLITZ matrix : 



% = 



to t\ • ■ • td-\ ^ 



\td-i • • • h t J 



we want to compute the signature of the associated Hermitian form. We didn't 
find any reference in the literature to this simple problem, although there are 
several methods proposed in the case of real Hankel matrices (see [12] and 
[32]). 

The signature of 7^ can be computed from the sequence of signs of its principal 
minors. The rule given by lOHVIDOV [18] and, independently, by one of us [33], 
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works even when some of these minors vanish. Once the minors are computed, 
the signature is obtained in 0(d) arithmetic operations. 



The problem is then reduced to the computation of the sequence of principal 
minors of T d . This can be achieved by the computation of the constant terms 
of the sequence of the symmetric subresultants of two polynomials as the next 
proposition shows. 

Proposition 12 Let T d be a Hermitian Toeplitz matrix, defined as above, 
and T the polynomial : 



T=-t- t ± x 



i d _ 1 x d - 1 + t d _ 1 x a + --- + t 1 x 



2d- 2 _ £j^2d-l 



with t 7^ and t = t + t. 

Let (5j)-i<7<2(2 be the sequence of symmetric subresultants of X 2d ~ l + 1 and 
T. For j = l,...,d, we have : 



where 5j is the j-th principal minor of T d . 

Proof : We can use Proposition 11 in this special case. But the result can 
also be seen directly as well. Indeed, we have for each j — 1, . . . , d : 



S,-(0) 



1 1 

i ; 

o 1 

-t ■ ■ ■ —tj-2 t 

-i : : 



-f t 



3-1 



j-1 
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1 o 



1 ••• 

—t —t\ • • • —tj-1 to tl • ■ • tj_l 

; /, ••• ; 

—t tj-i ■ ■ ■ t\ to 

□ 

Using FSSR Algorithm, we can then compute the signature of a Hermitian 
Toeplitz matrix of order d in 0(d\og(d) 2 loglog(d)) arithmetical operations. 

Brunie in [4] has shown that it is possible to improve the algorithm also to 
get the rank of the matrix, but this extra computation has an arithmetical cost 
of 0(d 2 ) operations. There still exists no fast solution to the rank problem. 



5.3 Toeplitz linear systems 

We now consider a much more popular application than the signature prob- 
lem. Let Td be a Toeplitz matrix of dimension d. Suppose it is invertible 
and we want to compute T^ 1 . Several authors have given fast algorithms to 
solve the problem. Brent, Gustavson and Yun in [2] have a solution us- 
ing PADE approximants, continued fractions and Euclidean algorithms. Their 
solution has a cost of 0(d\og(d) 2 loglog(rf)) arithmetical operations and uses 
the Gohberg-Semencul formulae. More recently Gemigniani in [13] and 
[14] has used the Schur decomposition of a matrix with the advantage that 
in defective cases no extra computation is needed. Both algorithms have the 
same cost. Bini and Pan give in [3] the state of the art on this problem. 

The solution developed here also works with the formulae of Gohberg- 
Semencul. However we use the symmetric subresultants; therefore we are 
able to manage the defective cases directly with the FSSR algorithm without 
extra computation. Our cost is the same as in [2] , although, in defective cases, 
we approximately divide computation time of by a factor two. Furthermore, 
our algorithm is fraction free, until the last step. 
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As it is one of our tools, we recall first the Gohberg-Semencul formulae 
[15]. 

Theorem 13 Let T d = (^-j)o<i,j<d-i be an invertible TOEPLITZ matrix. We 
denote by x = (xo, . . . ,Xd-iY the first column and by y = (yo, ... ,3/^-1)* the 
last column of T d l . If xq 7^ 0, we have : 



7-- 1 — J_ 
d ~x 



( 



x 



XQ J 



Ud-i 




yo ■ 



yo 



\ 



V • • • y d _i J 

... \ / 



• • • x x 



\Vd-2 • • ■ yo o) yo 







(*) 



Ifx = 0, there exists an extension T d+ i = (ti-j)o<i > j<d ofT d which is invertible 
and such that the first column of T d + X , say x = (x , . . . , x d ) , has its first 
coordinate different from zero. Let y = (y , . . . ,y d ) denote the last column of 
T d + X . In this case, we have : 



7— 1 — J_ 

d ~x 



t x 



\X d -l 



0) 




'yd 
















Xq J 




1° 


yo 








yi 



V • • • y d J 



( n ... n\ / 



: '•• 

\yd-i yo) 



x d 




Xi 



■ x d J 



(**) 



Therefore, \lT d — {U-j)o<i,j<d-i is an invertible Toeplitz matrix, the prob- 
lem is reduced to the computation of the vectors x and y or x and y depending 
on the situation. We can use the symmetric subresultants algorithm for this 
task. 
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Let us define the two polynomials : 




So = T^ 5 = -t--t-!X t_ d+1 X d - l + 1 X d +5X d+l +t d _ l X d+2 +- ■ -+t+X 



where coefficients t + and t_ are different from and satisfy t + + 1_ = t . The 
complex coefficients 7 and S will be determined later on during the computa- 
tion in order to apply Theorem 13. 

One can note that from F = S-i and G = S we can rebuild the matrix T 
using Proposition 11 : we have T = T d (S- 1: S ). 

Let (Sj)-i<j<2d+i be the sequence of symmetric subresultants computed with 
S-i and So- As T d is invertible, we have Sd(0) = (— l) d det(7^) 7^ (use Propo- 
sition 11). We will write S d = J2f=o SiX 1 . There also exist two polynomials 
U d -i = J2to u t X l and V d ^ = Ef=o ^X\ such that : 



X d - l S d = U d ^(X 2d+1 + 1) + V d ^T 7 , s . 



This relation can be translated into matricial terms as follows : 



/ 



1 







\ 



'• 







••• 1 











+ 











\Ud-l ) 



1 ••• 



••• ■•• : 



; ••. ••. 



Vo... 1J 
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-t-1 







—t-d+1 


■ ■ ■ 


-t- 


-7 
5 


—t-d+1 ' ' ' 


-t-1 


ti 


td-1 


-7 
5 


t+ 




td-i 



V 



o 



o * 



+ 7 



( o N 



o 

so 

Si 



\Vd-l) 



Sd+1 



V / 



d 



} d + 1 . 



If we subtract the first d lines from the last d ones, we obtain : 









with so = S'd(O) = (— l) a! det(7d) ^ 0. Therefore, we see that 



s 



is the 



\ v J 

first column of T d l . The same trick applied to T d gives the last column of our 
matrix. If v^-i ^ 0, we can apply the first formula of Gohberg-Semencul 
to conclude. 

By the proof of Lemma 1, we get v,i-i = —Sd-i(0) = (— ^) d det(7^-i). If 
Vd-i = 0, we have to compute the next symmetric subresultants, S d+ ±. There 
exist two polynomials, U d and V d , of degree at most d, such that : 



X d S d+1 = U d {X 2d+1 + 1) + V d T, 



7,<5- 
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In this case, deg(V d ) = d, because co d (V d ) — v d — (— l) d+l det(T^) ^ 0. If 
5d+i(0) 7^ 0, we see, by the same computation as in the generic case just 
above, that the coefficients of — V d / S d+ i(0) determine the first column of the 
inverse of : 



d+1 





7 


T d 


t-d+i 




t-i 


5 t d -i • • • ti 


to 



Therefore we have to choose the coefficients 7 and S in order to satisfy S d +i (0) = 
(-l) d det(T, +1 )^0. 

Proposition 14 Using the above definitions, suppose that det(7^_ 1 ) = and 
det(T d ) 7^ 0. Define the three vectors of dimension d : 



V_ 



1 ^ 

t-d+1 



v. 



< ^ 



V *1 / 



and e 







Then, the determinant of T d+ i satisfies : 



det(T d+1 ) = - det(T d ) • frV+'T^e,, + Se^V. + V+X^V- - t ). 

Furthermore, in the above relation, the coefficients \ + t T d ~ 1 e and e t T d ~ 1 \ '_, 
of '7 and 5 respectively, cannot vanish. 

Proof : We can factorize T d+1 as follows : 





o\ 




\ 






Id 


r 











\ $ td-i ■ • • h 


// 


v ••• 


l ) 
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with r = T d l 



( 7 \ 



-d+l 



= T d - 1 ( 7 e + V_) and : 



/ = to-(5e + V + ) t -V 1 ( 7 e + V_) 



Then, we have : 

/ = t - ( 7< 5 • eo'T^eo + 7 • V+'^eo + 5 • eo'T^V- + V + %" 1 V_). 



But £o l T d 1 e Q is, up to the factor 1/ det(T^), equal to det(7^_!) which is zero. 
Therefore, we obtain the stated formula. 



We know that T d is invertible; let (xo, . . . , Xd-if be the first column of its in- 
verse. Since det(7^_i) = 0, we have x = 0. If we suppose that 



V + % ^0 = 0, we have (0, t d . u 



( N 



Xi 



\Xd-l J 



= 0, and we can write : 










t-d+i 




t-i 


t d ^ • • • h 


to 



I N 



Xi 



Xd-1 



V u 7 











to 


t-i • 


■ • t-d+i 


tl 






td-1 














\ I o \ 



V / 
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We therefore conclude that : 



X d -1 

, , 



0. 



However, as T d is invertible, the equation T d • X = has only one solution, 
that is the zero vector. This leads to a contradiction since x±, . . . , a^-i are not 
all equal to zero. Therefore, the coefficient V + *7^ _1 e cannot vanish. A similar 
argument works with e *^ 1 V_. □ 

Now we are able to choose a pair (7,5) such that det{T d+ i) ^ 0. In fact, as 
the set of pairs (7,5) that make det(7^ + i) zero is a line, after three attempts 
we are guaranteed to find an acceptable value (for example, we try (0, 0), then 
(0, 1) and if, with both values, the determinant is zero, we can then use (1, 0) 
as a good coefficient). 

Before we describe the algorithm for fast inversion of a Toeplitz matrix, we 
have to make some important remarks. 

First, the polynomials L^-i and V^-i defined by: 



X d - l S d = U d ^(X 2d - l + l) + V d ^T^ (J) 

are obtained from FSSR applied to X 2d+1 + 1 and T 7i<5 with r = d + 2. As 
S d (0) 7^ 0, if degS d = d + 1, there exists k e such that k e = d. We can then 
compute M d . The coefficients on the second line of this matrix, M d , are exactly 
U d -\ and Vd-i, as we can see from the proof of Proposition 5. 

Otherwise, if deg S d < d+ 1, we observe that for the biggest £ such that ki < d 
we have the pair (Sk t , S^ e + 1) right-defective (indeed Theorem 3 shows that 
all other situations lead to Sfc(0) = for kp < k < kp + \). We know that in 
this case Sk t + 1 and S d are proportional ; the coefficient of proportionality is 
given by Theorem 3. From FSSR we obtain only : 



X ke Sk e +i = U ke (X 2d 1 + 1) + V ke Ty : s, 
Multiplication by the right coefficient provides formula ($). 

Furthermore, whatever the situation might be, in this call to FSSR, 7 and 8 
do not occur because we use a truncation to the order d — 1. 
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This provides the first column of T d 1 . The same computation applied to 
X 2d ~ 1 + 1 and Sq gives the last column. 

Next, we do not need any extra call to FSSR when we test, for example, 
(7,5) = (0,0), (1,0) or (0,1). The computations are different only for the 
last step, the transition from S d to S d+ ±, and we do not need to begin again 
the computation from S-i and So. This is the first advantage of our FITM 
algorithm over the one in [2]. A second advantage is that it is fraction- free. 

Finally we can rewrite our result in a TOEPLITZ-Bezoutian form. If (U, V) is 
a pair of polynomials of degree at most d such that 

X d S d+1 (S^, S ) = (X 2d+1 + l)V + UP, 
and if (u, v) is a pair of polynomials of degree at most d such that 

X d S d+1 (S_ 1 , S*) = (X 2d+1 + l)v + uP, 
then, in the non-degenerative situation, we have : 

Bez(U*,u)T d (S- U So) = S d (0)S d+1 (0)I d 

where I d is the identity matrix of order d. (It comes from a well-known matrix 
representation of Bezoutian - see [3], p. 156.) 

There are certainly relations between our computations and those proposed 
by Gemigniani in [13] and [14]. Bezoutians are used instead of symmetric 
sub- resultants. But, these algorithms start with quite the same polynomials. 
In the literature one finds several links between resultants and Bezoutians (see 
for example [20]). However, in our particular case, the relation between these 
two methods is not easy to describe and will be the object of future work. 

Of course, all that we have said in this sub-section can be simplified in the 
case of a Hermitian Toeplitz matrix. It has been described in detail in [4]. 

We can now summarize our results in the FITM algorithm for fast inversion 
of a Toeplitz matrix. 



6 Conclusion 

We have generalized the concepts introduced for the improvement of the 
Schur-Cohn algorithm. The sequence of sub-resultants defined for a pair 
(P, P*) can now be computed for a general pair of polynomials and the fast 
algorithm designed in the previous situation has been extended. 
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ALGORITHM FITM 


1JN rUl : 


Td = (ti-j)o<ij<d—i, a Toeplitz matrix of dimension d 


OUTPUT : 


T^ 1 if Td is invertible and, if not, a message that Td is not 




invertible 


TTVTTrf-lT A 1 TP ATIM\ 


C v2d+l i 1 




C„ rp j- j. v j., Y~d— 1 i j. i 

-->0 — J 0,0 — I t_iA I_^ +1 A ttd-l A + 




■ ■■ + tX 2d+ \ 


MAIN PART : 


- FSSR(5_i,5 ,d + 2) 




% we get Mfcj with 




% k\ the largest index such that k\ < d. 




- if ki = d and 5^(0) = or if k\ < d, and S^+iiO) = 0, 




Td is n °t invertible. STOP 




- compute Vd-i from and possible use of Theorem 3 




- FSSR(S_i,S i *,(i + 2) 




% we get Sr ,U k i, with 




% fc; the largest index such that k\ < d. 




— compute Vd-i from and possible use of Theorem 3 




Tf f \ d(V \l i -, — A 1 i~ h on ' / ^ ic nnmnntofl in q fnirnn 1 q ^"^^ 

11 vicg, 1 — (X 1, Llieil 1 j lo CUillJJULeLl Vldj lUIIHUld ^ J 




- If deg Vd-i = d — 1, then is computed via formula 




I ) 




— ii oeg v d— l \ o — J- ana aeg Vd—i <« o — compute o^-i-i 




using Mfc ; . 




— ii o^-i-i^uj /= u, rnen 1^ is compucea via iormuia ^ j 




— otherwise redo the computation of Sd+i with To i or with 




Ti, . 




% one of them will give 5^+1 (0) / 0. 


END. 





The effectiveness of the algorithms presented has been studied in [4] where 
they have been effectively programmed in TP language, using the DFT. It 
has been shown that the bounds are effective and that, for polynomials of 
degrees greater than 300 and coefficients bounded by 2 32 , these algorithms are 
faster than their counterpart programmed without DFT. 



Of course, the fast version of the Schur-Cohn algorithm has not changed, 
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but we can present applications to Toeplitz matrices which are new. It would 
be an interesting study to compare the different algorithms for the inversion 
of Toeplitz matrices and to explore the links between them. 
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